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Abstract 

Singlet pairing in a Fermi superfluid is frustrated when the amounts of fermions of each pairing 
partner are unequal. The resulting 'imbalanced superfluid' has been realized experimentally for 
ultracold atomic gases with s-wave interactions. Inspired by high-temperature superconductivity, 
we investigate the case of d-wave interactions, and find marked differences from the s-wave super- 
fluid. Whereas s-wave imbalanced Fermi gases tend to phase separate in real space, in a balanced 
condensate and an imbalanced normal halo, we show that the d-wave gas can phase separate in 
reciprocal space so that imbalance and superfluidity can coexist spatially. We show that the mech- 
anism explaining this property is the creation of polarized excitations in the nodes of the gap. The 
Sarma mechanism, present only at nonzero temperatures for the s-wave case, is still applicable in 
the temperature zero limit for the d-wave case. As a result, the d-wave BCS superfluid is more 
robust with respect to imbalance, and a region of the phase diagram can be identified where the 
s-wave BCS superfluidity is suppressed whereas the d-wave superfluidity is not. When these results 
are extended into the BEC limit of strongly bound molecules, the symmetry of the order parameter 
matters less. The effects of fluctuations beyond mean field is taken into account in the calculation of 
the structure factor and the critical temperature. The poles of the structure factor (corresponding 
to bound molecular states) are less damped in the d-wave case as compared to s-wave. On the BCS 
side of the unitarity limit, the critical temperature T c follows the temperature T* corresponding 
to the pair binding energy and as such will also be more robust against imbalance. Possible routes 
for the experimental observation of the d-wave superfluidity have been discussed. 



I. INTRODUCTION 



In a metal exhibiting superconductivity at low temperature, the amount of spin-up and 
spin-down electrons are equal, and the electron-phonon interaction leading to Cooper pairing 
has a given, fixed strength. Both the population of the spin-components and the electron- 
phonon interaction strength cannot be arbitrarily tuned, and this restricts the experimental 
study of superconductivity to some given values in parameter space. Nevertheless one would 
like to access a much larger region of parameter space to gain insight in pairing and the 
superconductivity. 

Superfluid Fermi gases have recently gained a lot of interest, precisely because of the 
accurate adaptability of the system parameters. The interaction strength between the two 
hyperfine spin states is an adjustable parameter. This allows to probe pairing and super- 
fluidity in the crossover between a Bardeen-Cooper-Schrieffer (BCS) state of weakly bound 
Cooper pairs and a Bose- Einstein condensate (BEC) of tightly bound molecules [l|, More- 
over, in a mixture of two hyperfine spin states of a fermionic element, the amount of each 
hyperfine spin component can be controlled experimentally. This permits to investigate 
the effect that a population imbalance between the spin components has on pairing [J H[ . 
Not surprisingly, these recent experimental breakthroughs [l|, 0, H, 0] have relaunched the 
theoretical efforts to understand imbalanced Fermi superfluids in the crossover regime 0]. 

The first theoretical study of Cooper pairing in an imbalanced Fermi mixture was per- 
formed in the context of BCS superconductors by Clogston p , who showed that a population 
imbalance destroys the superconductivity when the imbalance in chemical potentials is of 
the same order as the 'balanced' order parameter. Experiments confirm that imbalance 
frustrates pairing, and reveal that the excess spin component is preferentially expelled from 
the superfluid [U, 0]: demixing occurs 0]. More exotic pairing scenarios have been pre- 
dicted, most notoriously the 'Fulde-Ferrell-Larkin-Ovchinnikov' scenario jl] in which the 
Fermi spheres of the two components spontaneously deform, leading to Cooper pairs with 
nonzero center-of-mass momentum. 

When the temperature is raised and excitations are populated, the superconductivity may 
be restored by creating a 'balanced' pair condensate with an 'imbalanced' gas of excitations. 
This may lead to 'reentrant superconductivity' as proposed by Sarma [9(. In the Sarma 
state, the excess spin component is expelled from the superfluid, not in position space, but 
in energy space. 

In the context of the Sarma state, the case of an imbalanced Fermi gas with a ci-wave 
order parameter is particularly interesting. In the current experiments on superfluid Fermi 
gases, the temperatures are low enough so that only the s-wave partial wave matters, and 
the ci-wave scattering is much weaker than the s-wave interactions. However, the ci-wave 
order parameter has directions in momentum space where it vanishes, even at zero tempera- 
ture. This allows for a Sarma scenario where the excess spin component is expelled from the 
superfluid, not in position or energy space, but in momentum space. In this contribution, 
we show that ci-wave symmetry enables the superfluid to cope with imbalance all the way to 
temperature zero, using a similar scenario as proposed by Sarma for nonzero temperature. 
This leads to the conclusion that imbalance can stabilize the ci-wave pairs with respect to the 
s-wave pairs since the ci-wave superfluid is more robust against population imbalance. More- 
over, since also the ci-wave scattering length can be tuned through the Feshbach mechanism, 
we investigate the ci-wave superfluidity both in BEC and BCS regimes. 

The case of the ci-wave pairing in the BEC/BCS crossover is also interesting from the point 
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of view of high-temperature superconductivity [10J], where the order parameter is found to ex- 
hibit rf-wave symmetry (IH|. The current results, derived in the context of cold atomic gases, 
can also shed light on properties of the pseudogap state in the underdoped regime. This 



pseudogap (which appears to have the same ci-wave symmetry as the order parameter [121]) 
has been associated either with some competing order parameter in the normal state, or 
with the existence of "pre-formed pairs" . (l3| where a notable candidate is the non-condensed 
bipolaron[3]. Also in the current treatment, pre-formed d-wave pairs appear, and we present 
results for the pair binding energy of these objects as a function of temperature, density and 
interatomic interaction strength. 

Our formalism of choice to treat the imbalanced Fermi gases is path-integration. The 
path-integral formalism was effectively applied to study fermionic superfluidity in dilute 
gases using the approach with the Hubbard-Stratonovic transformation, choosing a saddle 
point, and performing the integration over the fermionic variables. This leaves an effective 
action depending on the saddle point value and the chemical potential. The effective action 
can be applied to study the Fermi superfluid in optical lattices [15| or to investigate vortices in 
Fermi superfluids [l6| . In Ref . [T7J , that approach was extended in order to take into account 
the fluctuations around the saddle point. The treatment starts from the partition function, 
which is the path integral over fermionic (Grassmann) variables. After the introducing the 
auxiliary bosonic variables and integrating over the fermionic variables, the exact expression 



for the partition function from [17| is the path integral over only the boson fields with 
an effective bosonic action. That action is then represented as a sum of the saddle-point 
contribution (which is calculated exactly) and the contribution due to Gaussian fluctuations, 
which is taken into account as a perturbation. At T = T c , this approximation for the 



fluctuations is equivalent to that of Nozieres and Schmitt-Rink [18 

The further development of this idea can be found, e. g., in Refs. 19j, |2l|. In Ref. 19 
the superfluid density is derived for a uniform two-component Fermi gas in the BCS-BEC 
crossover regime in the presence of an imposed superfluid flow, taking into account pairing 



fluctuations in a Gaussian approximation following Ref. [18j. In Refs. [201.1211]. the effects of 
quantum fluctuations about the saddle-point solution of the BCS-BEC crossover at T = in 
a dilute Fermi gas are included at the Gaussian level using the functional integral method. 
In Ref. [221, the superfluid density and the condensate fraction are investigated for a fermion 
gas in the BCS-BEC crossover regime at finite temperatures. The fluctuation effects on these 
quantities are included within a Gaussian approximation. The gas of interacting fermions 
in Refs. 19, 21, 22| is considered for the s-wave pairing and with no population imbalance. 



A study of the balanced cZ-wave system using the path-integral method can be found in 
Ref. [23j. Finally, the current work applies the path-integral theory to the imbalanced d- 
wave superfluids, including finite-temperature fluctuations, to show that <i-wave pairing is 
particularly robust against imbalance fluctuations. 

The formalism is presented in Section [III In Section IIII} we develop the mean-field ap- 
proach and discuss the resulting pair binding energy. In Section IIV| we also include the 
fluctuations to treat the finite-temperature case and determine the critical temperature for 
superfluidity. Near the unitarity limit where mean-field is known to fail it is necessary to 
include fluctuations, but also to incorporate the normal-state interactions in the description. 
The current approach achieves this through an expansion of the action around the saddle 
point that keeps terms related to the particle-hole excitations. To investigate these excita- 
tions, we calculate in Section [V] the structure factor for the d-wave and compare it to the 
structure factor in the s-wave pairing state. 
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II. FORMALISM 



A. d-wave interactions 

As in Refs. 0,|23, we start by writing down the partition function of the interacting 
Fermi g cLS cLS cL path integral over Grassmann variables: 

Z oc I V^n^V^n^expi-S) (!) 

Rather than use position and imaginary time variable, we have Grassmann fields tkk, n ,cr^'k,n,a 
that depend on the wave number k and the fermionic Matsubara frequency u n = ri7t//3 with 
n an odd integer and (3 = 1/(A;bT) the inverse temperature. Two different hyperfme spin 
states are trapped so that we include a spin quantum number a in the description. We'll 
denote the two states as 'spin-up', o =f and 'spin down', a =|. 

The action functional S = S + Si consists of a 'non-interacting part', Sq, and the 
interaction terms, Sj. The former is given by 

So = (~iuJn + k 2 - fJLtr) tp^ n ,a1pk,n,a (2) 

k,n a 

where \i a is the chemical potential fixing the amount of atoms of species a, and where the 
summations run over all possible indices of the Grassmann variables. We use units such that 
h = kp = 2m f = 1, where m/ is the mass of the fermionic atoms, and kp is the Fermi wave 
vector of the non-interacting, balanced Fermi gas with the same total number of particles. 
In what follows, we will use the average chemical potential fi = (/ij + /i^)/2 along with the 
difference in chemical potentials ( = — /i|)/2, rather than the chemical potentials of the 
individual species. 

The interaction terms of the action functional are written in a form that emphasizes the 
pairs of colliding fermions: 

S I = Yl Yl V PP ( k ' k ') + k >? +n,T^f-k,f -raV'f-k'.f -n'.jV'a+k'.f +n',T ( 3 ) 

q,m k,n k',n' 

Here V pp is the interaction potential. The wave numbers in this collision term are written 
as the sum of a center-of-mass wave number q and the relative wave numbers k, k' before 
and after collision. Similarly, the Matsubara frequencies are decomposed in a center-of- 
mass bosonic frequency Q m = 2mn/l3 and relative fermionic frequencies u n ,u n i. Here, we 
consider only interactions that couple fermions from different hyperfme spin states. We'll 



need a further assumption on the interaction potential to proceed. As in [23l . |24 
that the interatomic interaction potential can factorized as 



we assume 



V pp (k, k') = gT (k) r (k') (4) 
This is possible for s-wave pairing 

g = g s , T s (k) = 1, (5) 

and also for <i-wave pairing 



x, /i \ (k/hy /287T,, ,„ \ 



4 



Here, l2,o ($> V 9 ) is the spherical harmonic, and fei, ko are parameters fixing the range of the 
potential. The constant g < (g > 0) corresponds to attraction (repulsion). These constants 



can be related to the s-and c?-wave scattering lengths 25]. The usefulness of the factorization 



lies in the fact that it allows to rewrite the interaction terms as 

Si = g ^ ] ^4q,m^4q,m (7) 



q,m 

where we introduced the collective fields 



. _ ^ r (k) 

^q,m - "7^^f-k,f-na^|+k,f +n,T' 

n,k VP 

r (k) 

n.k 

Here is the system volume. The Hubbard- Stratonovic transformation can transform the 
product over these collective fields into a sum over them, at the expense of introducing an 
additional functional integration: 

z^Jv^v^Jv^v^i-S) (9) 



with the action 



S = 22 (~ iuJ n + k 2 - Ha) ^k,n,<7^k,n, CT (10) 



k,n,(r 



E/ A q , m A q , m -r ^ \ 

m,q V 5 1 / 

Note that the auxiliary fields Aq m ,Aq im are bosonic in nature, and characterized by the 
center-of-mass pair wave number and bosonic Matsubara frequency Q m = 2mTr/(3. The 
decoupling of the collective fields is necessary to perform the functional integral over Grass- 
mann variables, resulting in 

Z oc J T?Aq im X>Aq im exp I tr In [G _1 (q, to; k, n)] + - ^ A qim A qim j , (12) 

V " m,q / 

where — G _1 is the inverse Nambu tensor and the trace has to be taken over the fermionic 
degrees of freedom. 

The value where the (exponential) integrand becomes largest is called the saddle point. 
Interpreting A qm as the field of bosonic pairs, we can claim that when these pairs are 
condensed, the largest contribution derives from the terms with A 0j o = A. Performing the 
Bogoliubov shift, we change integration variables from A qm to 7 q>m where 



A q , m = Vy/3A5 m>0 5q,0 + 7q,m, (13) 
A q , m = V / Vy?A5 m)0 fl q ,0 + 7q,m- (14) 
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If, at this point, we choose the saddle point not as the q = state, but as a state with a finite 
wave number, equal to the difference between the Fermi wave numbers of each component, 
we obtain the FFLO state Since this has not yet been observed, we restrict the current 
calculations to q = pairs. This allows to split up the inverse Nambu tensor 

- G -1 (q, m; k, n) = <5 m , <5 q ,o (k, n)] + F (q, m; k) , (15) 

with the saddle-point contribution is 



"sp 



and the fluctuation contribution is 

We are left with the functional integration over the bosonic fields 7 qjm and jq >m . The simplest 
approximation consists in ignoring these fluctuations and setting G = G sp - this yields the 
saddle point results and will be explored in the next subsection. IIIII B. Expanding In [G _1 ] in 
successive orders of F yields a perturbation series in 7 qjm corresponding to an ever increasing 
diagrammatic expansion, with possible Dyson resummations. The term of order F 2 is still 
quadratic and we calculate it in subsection IIIIIC. Up to second order: 

Z oc exp {— S sp } x J 2}'7q,m^7q,m exp {— Sfi} (18) 

with 

S sp = tr In [&~p] - — AA (19) 



and 

S n = - tr (GJFGJF) - 
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Sfi = - tr (G sp FG sp F) ^ 7 q ,m7 q , m - (20) 



q,m 



Since the partition sum is a product Z = Z sp x Zfi , the corresponding thermodynamic 
potential will be a sum of a saddle-point contribution and fluctuations: F = F sp + Ffi. 
The contributions are defined by Z = e~@ F and Z sp ji = e~^ FspJl . These thermodynamic 
potentials will be necessary to calculate the two number equations in subsection IIIIID. 

B. Saddle-point action 

The trace over G" 1 (k, n) can be performed, yielding 

Z sp = exp (-S sp ) = exp {-(3F sp ) 
with S sp the saddle-point action 

S sp = -J2 ln - C - E k ) {-iu n + C - E k )} - — AA. (21) 

k,n ^ 
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where the reader is reminded that ( = — fi-\)/2 is the difference in chemical potentials. 
The Bogoliubov energy is 

(22) 



(k 2 -^y + \r(k)A\\ 



The sum over fermionic Matsubara frequencies can be calculated. We find for F sp , the 
saddle-point thermodynamic potential per unit volume, 



± sp 



dk 

(2^y 



(3 



In (2 cosh p( + 2 cosh (3E k ) - £ k 



i|A| 2 
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(23) 



where the fermion energy is £k = k 2 — fi. Thus, the saddle-point result is generic for all 
interaction potentials of the form 



C. Quadratic fluctuations 

When the terms of order O (F 3 ) and higher are neglected in Sfi, the functional integral 
over 7 q , m , 7q lTO in expression (|T5|) can be performed. The result is written as 



S 



fi 



(3V 



1 
2 



/ J^ ln [l-A^i^Cq^Q^)! 2 - |7Wi )2 (q,iQ^)| 2 ] (24) 



where now the trace is to be taken over bosonic Matsubara frequencies and center of mass 
wave numbers. Here, 



A* 1,1 (q,ifi n ) 



x 



dk 



(27T)' 



T 2 (k) 



2k 2 



+ 



sinh /?i£ k _s 



2£ k _ a cosh /3.E k _a + cosh f3( 



(iQ n - E k _a + £ k+ a J (-E k -f + £ k _a 



^ifi n + -Ek-f - -^k+f J (iQ n + ^k+f + ^k-f 



A (a) , (25) 



with the parameters A (a) which describe the coupling strength for the s-wave and <i-wave 
pairings [251] : 

1 2 
( a s) = z > {ad) = — r , (26) 



87ra s 



7ra 
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and 



M 1>2 (q, m n ) = - 1 A| 2 J ^r 2 (k) r (k + 1) r (k - 1) 

sinh f3E k _* 

X 3 

2£ k _s fcosh/3£ k _a +cosh/?C 



x 



1 



iQ n + -E k -a - -^k+a ) ( iO n + -^k-s + -E^k+a 



(27) 



In the particular case of the s-wave pairing and of the balanced fermion gas, (|25|) and 
( |27|) are equivalent to the matrix elements derived in Ref. j20|. In the present treatment, as 
distinct from Refs. (20I . 21]. we do not assume the low-temperature limit. 



D. Gap and number equations 

The gap equation is determined by S sp alone, through 5S sp /5A = 0. The gap equation 
can be written in a unified form for the s-wave and ci-wave pairings, 

/ ' r < k » 2 ( 2gk (J^ C0S MC) - if) + * « - ° W 

The number equations are determined from the thermodynamic potential through 

\dnJ TV> 

(-) 

\SUt,v, 

where n = + is the total local density, and <5n = n-j- — n-j- is the local population 
imbalance. For a finite temperature below T c , the chemical potentials ji and £, and the 
gap A are determined self-consistently as a solution of the gap equation ff28l coupled with 
the number equations (1291) and fl30|) . In principle, we can write the exact thermodynamic 
potential F = F sp + Ffi + F ot her where F sp and Ffi are given by expressions (1231 and fl24|) . and 
Fother comes from the contributions of all higher order terms, O (F 3 ) , in the exact action. 
The local density and the local population imbalance can be written as a sum of several 
contributions, 

n = n sp + n f i + n otherj (31) 
5n = 5n sp + 5n fl + 5n other , (32) 

where n sp and 5n sp are the saddle-point results, n// and Sn/i are the fluctuation contributions, 
and n ther, 8n ot h e r are higher-order fluctuation contributions to the density and population 
imbalance, which are neglected in the present treatment. 



-5n, 



(29) 
(30) 
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The saddle-point contributions to the density and to the population imbalance are ob- 
tained using the saddle-point term of the thermodynamic potential (123]) and Eqs. (1291) and 
(00}: 



1 



n 



sp 



5n 



sp 



2?r 2 
1 

2^ 



k*dk 1 1 - ^ 



sinh (fiE k ) 



k 2 dk 



E k cosh (PC) + cosh (pE k ) 
sinh(/?C) 



cosh (p() + cosh (pE k ) ' 



(33) 
(34) 



The fluctuation contribution to the number equations is determined on the basis of the 
fluctuation contribution to the thermodynamic potential: 



5n fl 



PJ (2-) 3 n : 

where the functions J (q, z) and K (q, z) for a complex argument z are given by 
1 



(35) 
(36) 



J (q, z) 
K(q,z) 



T(q,2) 

1 

r(q,z) 



aMi^ (q, g) aM 1|2 (q,z) 
AM*"*) ^ ^M*"*) ^ 

(q, -z) '-Ml^A _ ^ (q , _ z) ^1^) 



with 



(37) 
(38) 

T (q, z) = M hl (q, z) Mi,i (q, -z) - M h2 (q, z) M h2 (q, -z) . (39) 

The functions (q, 2) and Ai\ y2 (q, 2) of the complex argument z are analytical in the 
complex z-plane except the branching line, which lies at the real axis z = uj. Similarly to 
Refs. [l?], 26], the summations over the boson Matsubara frequencies in (|35|) and f[3"6"j) are 
converted to the contour integrals in the complex plane as described in the Appendix. Here, 
we write down the final result for the fluctuation contributions to n and 5n: 



n f i 
Sri fl 



dq 




(27T) 3 




dq 




(2vr) 3 





Im 



Im 



J(q,uJ + iff 

e /3(w+i 7 ) _ I 

K (q, uj + ry) 

P f3{u+i~f) _ 1 



1 

du + ~p 5^ ^ (q ; ^n) 



n=— no 
MO 



doJ + - ^ If (q, 



(40) 
(41) 



n=— no 



Here, the number n is chosen arbitrarily, and the parameter 7 lies in the range Q no < 7 < 

In particular, if one chooses no = 0, the formula (IA5I) leads to the expression for the 
fluctuation contribution to the fermion density similar to that derived in Ref. 171 ]: 



n f i 



dq 



(2tt)V 



dujS (q, cu) . 



(42) 
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Here, the structure factor is 



5 (q, w) = e ^ _ x , <3 -> +0. (43) 



The results obtained in the present section extends the path-integral approach of Ref. [21 



to the case of the ci-wave pairing and of an imbalanced Fermi gas at arbitrary temperatures. 



In agreement with the proof made in Ref. 26], the function 



Q (q, uj) = lim {Im [J (q, uj + id)] } (44) 

5^+0 

at T = T c is equal to zero at uj = 0. Furthermore, Q (q, uj) changes its sign as uj passes 
through uj = 0. This is necessary to ensure that the relative contribution to the fluctuation 
density from excitations with given (q, a;) remains positive; this contribution is proportional 
to S (q, uj). 



III. ROBUSTNESS OF THE d-WAVE PAIR BINDING ENERGY 

First we look at the saddle-point results for temperature zero, in order to investigate the 
pair binding energy. In the limit of temperature zero (3 — > oo, the gap equation becomes 

_ 2 _ f dk / 6(£ k >Q 1 
" n{k F a d ) 5 ~ J W( 2E V 2F 

with O the logical Heaviside function. Simultaneously the two saddle-point number equa- 
tions (|33l . (l34|) become 

/^{i-e ( *>o£} (46, 
p^e (£t < o. (47) 

The 0(i?k < C) function cuts off all the wave numbers with energy less than (. These are 
shown in Fig. dJ Near {k x , k y } = {v^2, \f2}kp the gap vanishes, and excitations are always 
present. To have non-zero imbalance in an s-wave superfluid, ( has to be of the order of 
|A|. This is the Clogston limit, and superfluidity will break down. However, for the ci-wave 
superfluid all values of £ lead to imbalance, and small values of ( do not destroy superfluidity. 

Solving for a given imbalance and interaction strength ad the saddle-point number and 
gap equations (|45j) - (l4Tj) . we can derive the saddle-point value A. This value is necessary to 
compute the fluctuation effects, but it has an interpretation by itself, namely as the pair 
binding energy. A corresponding temperature T* = |A| jks can be associated with the pair 
binding energy. In the BCS limit, superfluidity is destroyed by breaking up Cooper pairs. 
Thus, the transition temperature is determined by the binding energy of the Cooper pairs 
and T c rs T*. However, in the BEC limit, superfluidity is destroyed not by breaking up the 
bosonic molecules, but by phase fluctuations, and typically T c < T*. The BEC limit, with 
its tightly bound molecules, is relatively insensitive to the addition of atoms of one of the 
spin species: the imbalanced system can be described as a mixture of fermionic atoms and 
bosonic molecules. The BCS limit, however, is very sensitive to imbalance. Since in the 




1 

37f2 

1 5n 
3tt 2 n 
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FIG. 1: Even a small difference between the chemical potential leads to population imbalance. 
This plot shows the contours of E^ = £ for different values of £ in the k x , k y plane, for a d-wave 
interaction characterized by ko = k\ = 10/cf. Within the regions E^ < £, spin-polarized Bogoliubov 
excitations are present that carry the excess spin component of the imbalanced gas. 



BCS limit T c is directly related to the pair binding energy |A|, we gain insight in the effect 
of imbalance on s- and d-wave superfmids through the saddle-point gap. 

The result for |A| is shown in Fig. [2] for s-wave (top panel) and <i-wave (bottom panel) 
pairing. There are some notable differences between s- and <i-wave results. Firstly, for the 
d-wave interaction the x-axis is a function of l/(kpad) 5 in stead of l/(kpa s ). This means 
that the d-wave scattering length should be much closer to resonance as compared to the 
s-wave case before superfluidity enters the resonant regime. The absolute scale still depends 
on ko, related to the range of the interaction potential. Also the scale of the y-axis in the 
graph (representing |A|) has this dependence on the details of the potential embodied in 

r(k). 

A second difference between s- and d-wave resonant pairing, is that for s-wave interactions 
we find that there is pairing for all values of a s < 0. For <i-wave interactions, it is no longer 
true that for any attractive potential there is pairing. There needs to be a fatal attraction 
before pairing occurs on the BCS side. The BEC side, however, is more or less the same 
for s- and d-wave. Deep in the BEC regime it indeed should not matter whether we have 
s-wave or <i-wave internal parameter. 

A third difference, is that the ci-wave order parameter on the BCS side is much more robust 
to imbalance than the s-wave order parameter. For all negative scattering lengths, there 
exists a critical imbalance that destroys superfluidity in the s-wave system. However, in the 
d-wave case, there is a range of negative scattering lengths for which the pairs remain bound 
up to the maximal imbalance. This confirms our intuition that the excess spin component can 
be nicely stowed away in the minima of the gap, near the kjk F = {l/y/2, l/y/2} point. At 
these points, the gap vanishes naturally and it does not cost any energy to make excitations 
or to store broken pairs. 
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2.0 




10 4 /(k F a d ) 5 

FIG. 2: (color online) The dependence of the pair binding energy on the interaction strength for 
the s-wave (top) and d-wave (bottom) interaction is influenced by the imbalance 5n/n = 0.0 ... 
0.9. On the BCS side, the imbalance destroys s-wave Cooper pairs for all values of a s < 0, but 
fails to break up the cf-wave Cooper pairs. 

The (i-wave order parameter shows similarity to the s-wave when some imbalance already 
present (the 5n/n = curve looks like the s-wave curve for nonzero imbalance), but it is 
much more robust to imbalance. One can imagine increasing imbalance in such a way that 
it suppresses the s-wave pairing channel and still allows the c?-wave pairing channel. 

IV. CRITICAL TEMPERATURE FOR THE d-WAVE PAIRING IN THE REGION 
OF THE BCS-BEC CROSSOVER 

At finite temperatures, both phase fluctuations and amplitude fluctuations in A are 
important. The amplitude fluctuations dominate the thermodynamics in the BCS regime, 
whereas the phase fluctuations dominate in the BEC regime. This will be borne out in more 
detail by a study of the structure factor, in Sec. IV. For a given temperature T, density 
n and density imbalance Sn, we can solve the gap and number equations numerically and 
determine A, /i, (. The critical temperature can be found as the temperature where A 
vanishes. In Fig. [3l we plot the critical temperature in the case of the (i-wave scattering as 
a function of the inverse scattering length l/^kpad). The saddle-point results for the pair 
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breaking temperature T* = |A| /ks are plotted with the solid black curves, and the values 
of T c calculated taking into account the fluctuations are plotted with red full dots. 




-0.8 -0.4 0.0 0.4 0.8 1.2 



-0.4 0.0 0.4 0.8 1.2 



0.0 



FIG. 3: (color online) Critical temperature for the fermion system with the d-wave scattering 
calculated taking into account the fluctuations around the saddle point as a function of the inverse 
scattering length (red full dots). The saddle-point critical temperature T* is plotted with the solid 
black curves. 



For all considered values of the parameters ko, k\ of the d-wave scattering potential, we 
can see three following regions of l/(kFQ>d), with different behavior of T c . 

(1) A region corresponding to the weak-coupling regime (at l/(kpad) < 0). In this regime, 
with increasing l/(/cpa<i), the critical T c starts from the value T c = at a certain value 
l/(kpad), and rapidly increases. This is consistent with the finding in the previous section, 
that a critical strength of the interatomic interaction is required before pair formation occurs. 

(2) The region of the "plateau" around the unitarity point l/(fcpa<i) — 0, where T c varies 
extremely slowly. 

(3) The region corresponding to the strong-coupling regime (at l/{k F a d ) > 0), where T c 
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tends to the finite value T c fa 0.218. This is the same value as obtained in Ref. [17| for the 
s-wave case. Indeed we expect in the deep BEC limit the details of the internal structure of 
the molecule to be of secondary importance. 

Compared to the case of the s-wave scattering [l7| , the dependence T c (a d ) for the <i-wave 
scattering has a broad plateau around the point l/{k F ad) = both for the saddle-point 
results and for those taking into account the fluctuations. This plateau is explained by the 
fact that the factor l/{k F ad) enters the gap and number equations through its fifth power, 
1/ {k F dd) b , which varies very slowly as compared to the case of the s-wave scattering in the 
unitarity region. Another difference with the s-wave case is that there is a critical value of 
ad so that both T* and T c turn to zero. This means a minimum strength of the attraction is 
necessary to achieve pairing in the d-wave, whereas in the s-wave case for all values of the 
(negative) scattering length one has pairing. 

As the BEC limit is approached, T* strongly increases, while T c tends to a constant 
value. Again the behavior in the BEC limit is similar to that of the s-wave case, as can 
be expected. In the BCS regime, T c « T* as anticipated in the previous subsection. For 
k = 10, 5 and 3, T c is a slightly increasing function of l/{k F a d ) at positive l/k F ctd and 
remains everywhere lower than T*. For ko = 1, however, we see that T c achieves a maximum 
at a negative value of l/(k F ad) and then decreases to the BEC limit. This behavior of T c 
shows that at k = 10, 5 and 3, the anti-crossing of BCS and BEC regimes occurs at a d > 0, 
and that with decreasing ko, the region of the anti-crossing of BCS and BEC regimes shifts 
to lower values of the inverse scattering length. 



V. STRUCTURE FACTOR FOR THE s-WAVE AND d-WAVE PAIRINGS 

In the case with 7 — > +0, and for a balanced Fermi gas, the contribution n/i given 
by (HOj) is expressed through the integral (j4"2l with the structure factor S (q, uj) given by 
(14"3"|) . The structure factor is of particular interest, because it represents the spectrum of the 
elementary excitations of the fermion gas below T c . Further on, we analyze the structure 
factor S (q, u) at T = T c and the excitation spectra for the cases of the s-wave and d-w&ve 
scattering, comparing to existing results for the s-wave case [28]. Whereas the poles of the 
single-particle Green's function can be associated with single-particle excitations, the poles 
of S (q, ui) are related in the present case to the two-particle bound state. Note that the 
next term in the fluctuation expansion, proportional to the fourth power of A, gives rise to 
a spectral function the poles of which are related to the collective excitations of these bound 
modes. 

In Fig. UJ we have plotted the excitation region of the fermion gas in (q, cj)-space for the 
d-wave scattering at different values of the inverse scattering length l/a d and at T = T c . 
The solid black curve denotes the the lower bound Uq (q) = q 2 /2 — for the continuum of 
free two-particle excitations, which is determined by the inequality u > ujq (q). The dashed 
red curve corresponds to the solution of the equation u = (q), where fib (q) is the pole of 
the structure factor S (q,u). In the case when Qb (q) < ojq (q), i. e., when the pole fib (q) 
lies outside the continuum of two-particle excitations, Qb (o) is given by Qb (q) = ^>b (q) — 2yU, 
where Ub (q) is the energy of the two-body bound state [lTj. In the strong-coupling limit, 
the energy of the two-body bound state tends to u>b (q) = —E b + q 2 /2, where E b is the pair 
binding energy, which in this limit and at T = T c tends to (— 2/i). For a sufficiently weak 
coupling, the pole fib (q) lies within the continuum, and therefore the two-body bound state 
is damped. 
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FIG. 4: (color online) Excitation region of a gas of interacting fermions in (q, w)-space (the case 
of the c?-wave scattering) at T = T c for ko = 10, k\ = 10, and cos# = 1/2. The shaded area 
shows the continuum of the two-particle excitations. The black curve denotes the lower-frequency 
bound of the damping area. The dashed red curve shows the points given by uj = 0^ (q), where 
^6 (q) = (q) — 2/i with the energy oJb (q) of the two-body bound state. 



In the case of d-wave scattering, for all considered values of q, the two-body bound states 
are damped at l/a d < [Fig. @] (a, b)] and non-damped at l/a d <J [Fig. @] (d)]. For 
l/a d = [Fig. H] (c)], u (q) and fib (q) practically coincide. This allows to interpret, in the 
case of (i-wave scattering, the value l/a d = as the boundary between the regimes of the 
BCS-pairing (for l/a d < 0) and the BEC-pairing (for l/a d 0). 

The 3D plots in Fig. [5]represent the structure factors for the s-wave and d-wave pairing at 
weak-coupling. Because the scattering potential for the <i-wave scattering is angle- dependent, 
the structure factor 5" (q, a;) depends on three variables: S (q, uj) = S (q, cos 9, uj). Here, we 
discuss the results for S (q, uj) averaged over the directions, 

S(q,cu) = - [ S(q,co)sm9d9. (48) 

In Fig. [5] (a), the structure factor for the case of s-wave pairing is shown for l/a s = — 1, 
which lies on the BCS side of the resonance. In the BCS regime, the poles corresponding 
to the two-body pair excitations are damped since they lie in the continuum area. The 
spectral weight of those poles in the overall structure factor is significant only at small 
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FIG. 5: (color online) 3D plot of the structure factor S (q,co) for the s-wave scattering (a) and for 
the d-wave scattering (b) at T = T c in the weak-coupling regime. The critical temperatures are 
given in units of the Fermi temperature Tp = h 2 kp/ (2mfcs). At the top, there are the contour 
density plots for S(q,u>). 



wave vectors. This results in a peak around q = in Fig. [5] (a). Furthermore, there is 
a distinctive extremum of the structure factor at the boundary of the continuum of two- 
particle excitations. In Fig. [5] (b) we switch from s-wave to rf-wave pairing, but remain 
within on BCS side of the resonance: this panel shows the structure factor for the (i-wave 
pairing is for = —0.17. This value of the (i-wave scattering length is close to the lowest 
value of the inverse scattering length at which pairing can occur. Also in the (i-wave BCS 
regime, the two-body bound states are damped and therefore the peak corresponding to the 
two-body bound excitations has a finite width. However, the width of that peak in the case 
of the (i-wave scattering is relatively low. From this we can see that in the case of the d- 
wave scattering, the two-body bound state plays a significant role even in the weak-coupling 
regime. For the c?-wave scattering, as distinct from the s-wave scattering, the BCS pairing 
mechanism can be realized only in a narrow range of the inverse scattering length close to 
the lowest value of 1/ad from those, for which pairing can occur. 

Fig. [6] describes the strong-coupling case (on the BEC side of the resonance) where there 
is a non-damped isolated pole in the structure factor. Therefore, the structure factor in the 
strong-coupling regime contains a 5-like peak, which lies outside the continuum of free-pair 
excitations. In order to visualize those 5-like peaks, we use a finite damping parameter 
7 = 0.01. In the strong-coupling regime, the regular part of 5* (q, uj) is negligibly small with 
respect to the main contribution due to the aforesaid isolated pole, which describes the BEC 
pairing. 

VI. DISCUSSION AND CONCLUSIONS 
A. The unitarity limit 

A first point to discuss, is the predicament of mean-field theory in the unitarity limit, 
1/a — > 0. It is crucial that fluctuations are taken into account, as we have done in the 
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FIG. 6: (color online) 3D plot of the structure factor S (q,uj) for the s-wave scattering (a) and for 
the d-wme scattering (6) at T = T c in the strong-coupling regime. 



previous section. Higher-order fluctuation contributions can be taken into account, and 
many different approaches were developed for the balanced s-wave Fermi superfluid. These 
approaches differ in the types of higher-order processes they take into account. Strinati and 
co-workers [2t| work diagrammatically to improve on the results obtained by Nozieres and 



Schmitt-Rink for T = T c . Levin and co-workers 30] construct a finite temperature theory 
similar to that of Strinati and co-workers but include different diagrams in the summation. 
Alternatively, quantum Monte Carlo simulations [31] can be used to obtain results on the 
crossover physics the balanced s-wave Fermi superfluid and put the crossover theories to the 
test . At the unitarity limit the existing theories do not succeed to find the free energy of 
the system with more than 10% accuracy with respect to the Monte Carlo results. As such, 
we expect the current theory to have a similar level of accuracy in the unitarity limit for the 
(i-wave system. 

The problem that lies at the root of the difficulty to make a theory for unitarity is 
that one needs to take into account not only the fluctuations of the order parameter but 

are 



also the normal-state interactions correctly. The diagrammatic approaches |29|, [30 



based on a zeroth-order decoupling that emphasizes pair formation rather than normal-state 



interactions. Put in the language of functional integration [2l|], we have made a particular 



choice for the Hubbard- Stratonovic decoupling: we took ^ and ipip types of products of 
Grassmann variables. This emphasizes pairing, but when the pairing goes to zero at the 
saddle-point, the resulting normal state has no interactions, and fluctuation corrections 
are needed to remedy this. We could have made the choice to group and ipip and 
apply the Hubbard-Stratonovic scheme to decouple the four-product in these densities rather 
than in the pairs. The resulting saddle-point approximation would yield the random phase 
approximation (RPA) results for the interacting normal state; but it would lack pairing. The 
inability to include -on the level of a saddle-point approximation- both pairing and RPA 
normal-state interactions through the introduction of two collective fields is discussed by 
Kleinert [32| . who proposes variational perturbation theory as a solution [33|. The fluctuation 
expansion used in the current work goes beyond that of Ref . [2l| , in that we take not only the 
particle-pair tptp and hole-pair tptp excitations, but also particle-hole terms ipip are present. 
These contributions appear in terms that do not vanish as the saddle point goes to zero, 
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A — > 0, so that the normal state in the present treatment is the interacting Fermi gas rather 
than the ideal Fermi gas. As such, the present treatment will be better suited in the unitary 
limit. 



B. Routes for experimental observation 

Magnetically tuning the population imbalance in a Fermi superfluid is out of reach at 
present in high-temperature superconductors. In cold atomic gases it has been successfully 
demonstrated and applied to reach the superfluid regime. The currently realized atomic 
Fermi superfluids have s-wave symmetry of the order parameter. The <i-wave coupling 
strength generally is too small to dominate the s-wave scattering at low temperature. This 
can be overcome using a Feshbach resonance in the <i-wave scattering channel. c?-wave 



Feshbach resonances have been observed, for example, in various isotopes of rubidium|35 
However, to reach the unitary limit for the (i-wave scattering one needs a better control over 
the magnetic field than in the s-wave case, since the interaction parameter A, expression 
f l26|) . scales clS CL^ dbS compared to a s 1 . The current results suggest a different route towards 
(i-wave superfluidity: imbalancing the gas. The dominant s-wave pairing is easily suppressed 
on the BCS side of the resonance by adding imbalance, whereas (i-wave superfluidity is less 
sensitive to imbalance. Both the use of a (i-wave Feshbach resonance to obtain a large enough 
(i-wave coupling strength, and of imbalance to suppress s-wave pairing, will be needed to 
realize (i-wave superfluidity in the atomic gases. 

In an inhomogneous trapping potential, phase separation can occur in real space. For s- 
wave superfluidity, this leads to a balanced superfluid at the center of the trap, surrounded 
by a halo of imbalanced (or fully polarized) normal gas j34[. In effect, the excess spin 
component has been expelled from the balanced s-wave superfluid. The additional energy 
cost of placing the excess atoms fhigher up the trapping potential is compensated by the 
energy gained by allowing the balanced superfluid state to form. In the (i-wave superfluid, 
this energy balance is different. Increasing imbalance in the BCS side does not strongly 
reduce the free energy of the superfluid. As can be seen from Fig. [2J the (i-wave order 
parameter is not strongly affected. Therefore, expelling the excess atoms to the edge of the 
trap raises the total energy and we do not expect real-space phase separation. The situation 
is different on the BEC side: here, the free energy of the (i-wave superfluid is reduced by 
imbalance, and it becomes energetically favourable to expell excess atoms. 



C. Exotic pairing scenarios 

Note that the action for the fluctuations depends (through M.\\ and M.2\) on the choice 
of saddle point A. This means that the spectrum of excitations (obtained from the diago- 
nalized fluctuation action) also depends on the choice of the saddle point. Excitations for 
a vortex condensate may be different from excitations on top of a ground-state condensate. 
At nonzero temperature those excitations will be populated through Bose statistics. But 
the physics is more complex than just Bose populating excitations: the excitation spec- 
trum itself (the dispersion and lifetime of those excitations) is temperature dependent: new 
single-particle and collective modes appear and shift as a function of temperature. At zero 
temperature the only single-particle excitations are E^, the energy spectrum for breaking 
a Cooper pair, but at finite temperature, we also have the excitations of the thermal gas. 
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These consist in taking the atoms of a broken Cooper pair, and giving those atoms an extra 
kick: E'k+q - -^k- Besides those single-particle excitations, we will have collective excitations 
whenever M11M22—M21M12 = 0. 

In the case of an imbalanced Fermi gas, an alternative choice for the saddle point is 
Ae lkr where k represent a shortest wave vector connecting the Fermi surface of the minority 
component to that of the majority component. The resulting equations describe the Fulde- 



Ferrell-Larkin-Ovchinnikov state |8|, [36j. However, this state has not yet been reported 



experimentally, so we have restricted the present analysis to the usual pairing scenario. 



D. Conclusions 

We have investigated the imbalanced d-wave Fermi gas, both at zero and at nonzero 
temperatures, and as a function of the (i-wave interaction strength. We find that in the BCS 
regime, the (i-wave pairing is more robust to the presence of population imbalance than 
the s-wave case. For a range of interaction strengths, we find that the s-wave superfluidity 
is suppressed whereas the (i-wave superfluidity is not. This is shown to be related to the 
possibility of creating a polarized gas of excitations in the nodes of the gap. Rather than 
phase separation in real space, phase separation can occur in reciprocal space. An additional 
difference with the s-wave BCS case, is that a critical attraction strength is needed in the 
(i-wave case before pairing can occur (in s-wave pairing occurs for all attractive interac- 
tion strengths). In the BEC regime, the symmetry of the pairing interaction plays a less 
important role: as the molecule gets more tightly bound, the details of its internal wave 
function matter less, and we retrieve known results for the s-wave system in the same BEC 



limit [28|, |29|, |30[. We then investigate how our results are affected by Gaussian fluctuations, 
important both to describe the nonzero-temperature thermodynamics. Both the critical 
temperature and the effect of temperature on the spectral density of the excitations are 
calculated. Our investigation of the structure factor reveals that for the (i-wave scattering, 
the damping of the pole for S (q, u) is very small even in the BCS regime, in contrast to 
that for the s-wave scattering. The critical temperature in the BCS regime reflects the pair 
binding energy. This implies that the critical temperature for the (i-wave superfluid in the 
BCS regime will also be more robust against population imbalance. 



Acknowledgments 

This work has been supported by the FWO-V Project Nos. G. 0356. 06, G. 0115. 06, 
G.0435.03, G.0306.00, the WOG Project No. WO.025.99N, and the NOI BOF UA 2004. 

APPENDIX A: MATSUBARA SUMMATION FOR THE DENSITY 

Let us consider the contour integral on the contour C shown in Fig. [7J 

1 ^ f W-dz, (Al) 



2ni J c eP* 

where the points z = iQ n with |n| > no lie inside the contour, and the other points z = ifl n 
are outside the contour. The function / (z) possesses the following properties: (i) it is 
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analytic in the entire complex z-plane except, possibly, the branching line on the real axis, (ii) 
/ (z) decreases at Re z — > — oo faster than so that the integral f° / (u ± ij) duj, where 
uj and 7 are real, converges. The functions J (q, 2) and K (q, 2) determined, respectively, 
by Eqs. ( 1571) and (J38J), satisfy these conditions. The fraction e ^\_ x has the poles at z — ifl n , 
n = 0, ±1, ±2, . . .. The residues of e)3z 1 _ 1 in the points z = iQ n are equal to i. 



Imz 





FIG. 7: (color online) Integration contour in the complex z-plane. The full dots indicate the poles 
z = i£l n , n = 0, ±1, ±2, .... 



On the one hand, the integral flAll) is equal to the sum of the residues of the function 
in the points z = iVL n inside the contour C: 



,,3: 



On the other hand, the integral / is 

j = -L 

2m 



I E / 
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where the parameter 7 satisfies the inequality 

^n < 7 < ^n +l- 

It follows from the equivalence of ( 1A2I) and ( 1A3I) that 
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7T 



Im 
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duj + / (ifl r 



(A2) 

(A3) 
(A4) 

(A5) 



According to the theorem (1A5|) . the fluctuation contributions to the density and to the 
population imbalance can be represented as 



n f i 
8n fl 




Im 



Im 



K (q, uj + i'f) 

e /3(u)+i-y) I 



^ "0 

du) + - J (q, zO r , 
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du} + ~p E ^ (q>^«) ) • 



(A6) 
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n=— no 
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As follows from the above analytical transformations of the integrals in the complex z-plane, 
the sum (1A5I) does not depend on the choice of the number no and (for a given no) on the 
value of 7 within the range given by flA4j) . 
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